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f^ • We are investigating numerically the non-linear behaviour of a space- 

periodic MHD system with ABC forcing. Most computations are per- 
^f\ ' formed for magnetic Reynolds numbers increasing from to 60 and 

K^ ■ a fixed kinematic Reynolds number, small enough for the trivial solu- 

^^ , tion with a zero magnetic field to be stable to velocity perturbations. 

(•~^ I At the critical magnetic Reynolds number for the onset of instability 

t^ ' of the trivial solution the dominant eigenvalue of the kinematic dy- 

namo problem is real. In agreement with the bifurcation theory new 
steady states with non-vanishing magnetic field appear in this bifurca- 
^ ■ tion. Subsequent bifurcations are investigated. A regime is detected, 

_ ^ . where chaotic variations of the magnetic field orientation (analogous to 

^ \ magnetic field reversals) are observed in the temporal evolution of the 

(— I ' system. 
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1 Introduction 

In numerical modelling of physical phenomena two approaches are feasible. 
On the one hand, one can develop a model as detailed, as possible. Then an 
obvious difficulty is encountered: usually simplifications are necessary, owing 
to complexity of natural phenomena and/or lack of precise data. Neverthe- 
less, if the model is accurate enough, it displays behaviour similar to that of 
the underlying physical system, when fixed parameter values equal or close 
to those in the original system are employed. On the other hand, one can 
consider a class of approximate models, refiecting only global features of the 
underlying physical system. Parameter values are varied in order to study how 
this affects the overall behaviour of the system. The second approach is ad- 
vantageous in that properties characterising the class of systems are revealed; 
computationally it is usually less demanding, than the first one. 

We follow the second approach to examine a succession of temporal regimes 
which can occur in magnetohydrodynamic (MHD) systems, when Reynolds 
numbers are increased. In particular, we are concerned with two questions: 
Which features of MHD systems are responsible for magnetic field reversals - 
an attribute of the Earth's magnetic field? Which sequences of bifurcations 
can bring an MHD system to such a regime? 

A comprehensive review of low-dimensional ODE dynamo models of mag- 
netic field reversals of the Earth, discussed in literature prior to 1994, can 
be found in Jacobs (1994), Chapter 5. The models were categorised into two 
distinct classes. In the models of the first class it was assumed that reversals 
were related to MHD instabilities, triggered by finite-amplitude perturbations 
of otherwise stable nonlinear dynamos. The perturbations were supposed to 
be controlled by random processes, and various physical mechanisms of non- 
MHD nature were proposed as their possible sources. In models of the second 
class the reversals were a feature of the system, with no external perturbations 
required to initiate a reversal. The model of solar magnetic field reversals by 
Zeldovich et al. (1983) also falls within this class. Hide et al. (1996) and 
Moroz et al. (1998) generalised the Rikitake's disk model. 

Armbruster et al. (2001) proposed a 7-dimensional ODE dynamo model 
obtained by center manifold reduction with the use of symmetry consider- 
ations. They considered convection in a spherical shell without rotation for 
such parameter values, that mode interaction took place and heteroclinic cycles 
emerged in the phase space, connecting the 7 equilibria present in their system. 
The heteroclinic cycles were stable due to the symmetries of the system. In- 
troduction of magnetic field could transform the hydrodynamic equilibria into 
MHD ones, with a non-vanishing magnetic field. The resultant heteroclinic 
cycles involved changes of orientation of the magnetic field, i.e., reversals of 
the magnetic dipole (non-dipole magnetic components were ignored in their 
truncated system). 

In the model of Melbourne et al. (2001) reversals were also linked with 
an underlying heteroclinic cycle. To obtain the cycle, they considered con- 



vection in a spherical shell with rotation and derived a normal form near 
a codimension-3 bifurcation, assuming that the flow sustaining the dynamo 
possessed a certain group of symmetries. The normal form, truncated up to 
terms of the third order, involved 9 parameters. When several higher-order 
small symmetry-breaking terms were added in the normal form, for some pa- 
rameter values the behaviour of the system resembled the Earth's magnetic 
field reversals. 

Reversals were observed in simulations of MHD systems with the geom- 
etry and boundary conditions corresponding to those of the Earth in the 
3-dimensional model (Glatzmaier and Roberts, 1995, 1996; Roberts and Glatz- 
maier, 2001), involving the Navier-Stokes equation with the Coriolis, buoyancy 
and Lorentz forces, the magnetic induction and the temperature equations. A 
2|-dimensional model with a limited resolution in the azimuthal direction, de- 
rived from these PDE's, also featured magnetic reversals (Sarson and Jones, 
1999; Sarson, 2000). Magnetic field generated in these models was predomi- 
nantly dipolar; the lengths of time intervals of constant polarity and a short 
duration of each reversal were consistent with those of the Earth. However, 
the computations were very demanding in CPU time, prohibiting numerical 
identification of complex bifurcations in such systems. 

In simulations of Glatzmaier and Roberts reversals were robust. Glatz- 
maier and Roberts (1996) used a modified set of equations (compared to the 
one studied by Glatzmaier and Roberts, 1995: compositional buoyancy was 
added and the anelastic approximation was used instead of the Boussinesq 
approximation), but this did not result in significant changes of the temporal 
behaviour of magnetic field. Despite the model was run very far from the 
regime of the geophysical parameter values (Jones, 2000), there was a strong 
similarity between the model output and the geodynamo. This is another in- 
dication that their MHD attractor is structurally stable (i.e., persists when 
equations, geometry, or parameter values are considerably varied). Roberts 
and Glatzmaier (2001) considered three different radii of the inner core, cor- 
responding to its present size and some values in the past and in the future 
(geodynamo was driven by the energy influx due to the latent heat release 
at the freezing surface of the inner core, slowly growing in time). Again, in 
all the three cases magnetic field of comparable strength was generated, and 
the magnetic field exhibited a quantitatively similar temporal behaviour, the 
system in the future displaying a stronger temporal variability. 

2 Statement of the problem 

Equations of magnetohydrodynamics are invariant under the symmetry, which 
preserves flows v and changes orientation of magnetic fields b: 

/,;(v,b)^(v,-b). (1) 



One can anticipate the following scenario of development of a temporal regime 
involving magnetic field reversals. Suppose for certain values of a control pa- 
rameter the system possesses a stable steady state, in which magnetic field 
vanishes. When the parameter is varied beyond the critical value, two steady 
states with a non-vanishing magnetic field emerge in a pitchfork bifurcation 
(generic to systems with the symmetry (1) ). In subsequent bifurcations the 
steady states bifurcate into more complex attractors, related by (1) and re- 
maining so far distinct. As the parameter is further varied, the attractors of 
the system join (e.g., due to collision of attractors, or via heteroclinic connec- 
tion) and a unique attractor emerges. A sample trajectory on the attractor 
jumps between former attractors (which degenerate into unstable invariant 
sets). Since the jumps involve changes in the magnetic field orientation, the 
resulting intermittency can thus be linked with the reversals. 

The goal of the present study is to check whether this scenario is feasible. 
As a test case we consider the Navier-Stokes equation 

dv 1 

-_ = V X (V X v) + -Av - b X (V X b) + f - Vp (2.a) 

ot ri 

and the magnetic induction equation 

^ = V X (v X b) + ^Ab (2.6) 

at Km 

under the solenoidality conditions 

V-v = 0, V-b = 0. (2.c) 

Here v is velocity of the flow, b - magnetic field, p - pressure, R and /?„ 
are kinematic and magnetic Reynolds numbers. (Since in the problem under 
consideration characteristic length and velocity magnitudes are of order one, 
the Reynolds numbers are defined as inverse viscosity and inverse magnetic 
diffusivity, respectively.) The fields are supposed to be 27r-periodic in space. 
The force 

f = uo/i? (3) 

is assumed, where uq is an ABC flow 

uq = (Asinx3 -|- Ccosa;2, Bsinxi + Acosx^, C sm.X2 + B cosxi). (4) 

For this force, 

V = Uq, b = (5) 

is a steady solution to (2)- (4) for all Reynolds numbers. 

This particular system is considered for the following reasons: First, space 
periodicity enables one to use pseudo-spectral methods (see Canuto et ai, 
1989; Boyd, 1989), which are computationally efficient. Second, linear stability 
of the system (2)-(5) to hydrodynamic (Galloway and Frisch, 1987; Podvigina 



and Pouquet, 1994; Podvigina, 1999) and magnetic (Arnold and Korkina, 
1983; Galloway and Friscli, 1986; Galanti et al, 1992; Childress and Gilbert, 
1995) perturbations has been investigated both numerically and analytically 
(nonlinear regimes for some parameter values were also explored by Podvigina 
and Pouquet, 1994; Podvigina, 1999; Galanti et al, 1992; Feudel et al, 1996; 
Brummell et al, 2001). These studies provide guidance, in what parameter 
ranges the targeted type of behaviour may be observed. In particular, growing 
magnetic modes exist for the flow (3) in broad intervals of constants A, B and 
C. Third, since the critical magnetic Reynolds number is typically of the order 
of 10, only moderate Rm need to be considered, and a relatively low resolution 
(32^ Fourier harmonics) suffices. 

Magnetic field growth rates were calculated by Galanti et al. (1992) for the 
fiow (4) with the coefficients satisfying B = C, A^+B^ + C^ = 3, < B/A < 1 
for various values of Rm- For R^ = 12 there are three windows in B/A of 
positive growth rates (see Fig. 2 ibid.). We have checked that in the window 
0.7 < B/A < 0.9 the dominant eigenvalue is real, and in the other two windows 
it is complex. Distinct steady states with a non-zero magnetic field appear in 
a bifurcation of the trivial steady state (5), only if the magnetic induction 
operator has a zero (hence real) eigenvalue at the bifurcation point. This has 
suggested to set 

A = l, B = C = 0.75; (6) 

these values of constants in (4) are assumed throughout. 

With our goal in mind, in the present study we focus at non-hydrodynamic 
attractors. The hydrodynamic global stability of Uq guarantees that magnetic 
field does not vanish in saturated regime, if v = Uq is a kinematic dynamo. We 
have verified that for i? = 4, v = uq is a unique attractor of the hydrodynamic 
system ( (2.a,c) with b = 0), and most computations are made for this value 
of R. Several runs are also performed for higher R. 

3 The group of symmetries 

The group of symmetries without inversion of time of an ABC fiow with B = C, 
which we denote by ^, is comprised of two independent symmetries 

TT Vr TT 

Si : Xi -^ — — — ^3, X2 -^ —— — X2, X3 -^ — — — Xi, 
TT TT vr 

S2 ■ a:i — > - + X3, X2 -^ -- + X2, xs ^ - - xi, 
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Se = S1S2 : xi ^ -- + X3, X2 -^ - - X2, X3 ^ - + Xi 

S7 = Sisl : Xi -^ —Xi, X2 ^' n — X2, 0:3 -^ vr + X3, 

and the identity transformation sg = e (Arnold, 1984; Dombre et ai, 1986; 
Podvigina and Pouquet, 1994). The group is isomorphic to the symmetry 
group of a square, D4. 

Any symmetry of the force (3), (4) which does not involve inversion of time 
is also a symmetry of the Navier-Stokes equation (2. a) with b = 0. The group 
of symmetries of the MHD system (2)-(4), (6), which we denote by Ti., is a 
direct product of Q and the 2-element group Z2, generated by h (1). It has 16 
elements, which are either Sj, or hsi. 

Problems in hydrodynamics and magnetohydrodynamics (see e.g. the re- 
view by Crawford and Knobloch, 1991) often involve symmetries, induced by 
the geometry of the region where fluid is contained. In particular, if convection 
in a plane layer is restricted to space-periodic flows with equal periods along 
two Cartesian axes parallel to boundaries, and conditions at the boundaries 
are different, the group of symmetries of the system is isomorphic to Q. If in 
addition magnetic fleld is present, the symmetry group is isomorphic to Ti. 
Thus, though the system under consideration (2)-(4), (6) is highly symmetric, 
it can not be regarded as an exotic one. 

Symmetries are useful for categorisation of attractors. Let A be an attrac- 
tor of a dynamical system, invariant under a symmetry g: g{A) = A. Two 
cases can be distinguished: either A is pointwise invariant, i.e. gix.) = x for all 
points X G ^, or it is invariant only as a set, with g{'x) ^ x for some x G ^. In 
what follows, only the symmetries for which an attractor is pointwise invariant 
are regarded as symmetries of the attractor. 

4 Attractors of the MHD system for i? = 4 

Results of computations are summarised in Table 1 (see also the bifurcation 
diagram Fig.l). The magnetic Reynolds number has been increased step 1, 
except in the interval 37 < Rm < 40, where a high density of bifurcations 
has required smaller steps. For every considered Rm an attractor has been 
obtained in a run with an initial condition, which is a perturbation of the 
trivial steady state (5) with the energy of 10~^ in each Fourier harmonics 
spherical shell. The range of existence of attractors has been determined by 
continuation in parameter: the runs are done with an initial condition, which 
is a point on the attractor for close values of Rm- 

For 15 < Rm < 58 the system under consideration possesses multiple 
attractors. For 16 < Rm < 58 attractors, which we have detected, belong to 
two independent families; within each family they are genetically related by se- 
quences of bifurcations. Attractors from different families have distinct groups 
of symmetries. The families can be distinguished also by their time-averaged 
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Table 1. Attractors, detected for the MHD system with the forcing (3), (4), (6) 
for R = 4 and < Rm < 60. The third column shows the number of elements 
of the symmetry group for which an attractor is pointwise invariant, and the 
fourth column - generators of the group. Indices in the labels of attractors 
(second column) have the following meaning. For periodic orbits emerging in 
the period doubling cascade the superscript M refers to an Mr-period orbit, 
r being the period of the foremost orbit. The first subscript labels different 
attractors of a given kind. The second subscript labels mutually symmetric 
attractors. 



R 


Attractors 


Number 
of sym- 
metries 


Generators 


Eu 


Em 


Rm<U 


steady state 5*0 


16 


Si, S2, h 


1.06 





15 < Rrn < 24 


steady states Si^i {i = 1,2) 


8 


hSi, S2 


1.06-0.92 


0-0.035 


25 <Rm< 37 


periodic orbits Pi^i {i = 1,2) 


4 


hsi, S3(= si) 


0.94 


0.03 


37.5 <Rm< 38.5 


steady states Si^i {i = 1,2) 


8 


hsi 


S2 


0.95 


0.025 


39<R^< 39.2 


periodic orbits P^i {i = 1,2) 


8 


hsi 


S2 


0.98 


0.02 


39.3 <Rm< 39.4 


periodic orbits Pg^,* (^ = 1)2) 


8 


hsi 


S2 


0.98 


0.02 


Rm = 39.45 


periodic orbits Pg^^ {i = 1,2) 


8 


hsi 


S2 


0.98 


0.02 


Rm = 39.5 


chaotic Ci^i = P^^ {i = 1, 2) 


8 


hsi 


S2 


0.98 


0.02 


Rm = 39.6 


periodic orbits Pg^j {i = 1,2) 
& periodic orbit P3 


8 
8 


hsi 
hsi 


S2 
S2 


0.97 
0.95 


0.025 
0.03 


39.7 < i?„ < 41 


periodic orbit P4 


8 


hsi 


S2 


0.93 


0.035 


i?„ = 42,43,44 


tori Ti,, (i = 1, 2) 


4 


hsi 


S3 


0.93 


0.04 


Rn. = ^5 


chaotic 6*2,1 (^ = 1)2) 


4 


hsi 


S3 


0.93 


0.04 


Rm = 46, 47 


chaotic Cs^i {i = 1,4) 


2 


hsi 


0.94 


0.03 


Rra = 48, 49 


chaotic C^^i (i = 1,4) 


2 


hsi 


0.92 


0.035 


Rm = 50, 51 


chaotic Cs^i (i = 1,2) 


2 


hsi 


0.91 


0.04 


52 < Rm < 58 


toriT2,, (^=1,4) 


2 


hsi 


0.9 


0.05 


lQ<Rm< 25 


periodic orbit Pg 


4 


hS2 


0.4-0.7 


0.07-0.2 


26<Rm< 60 


torus T3 


4 


hS2 


0.3-0.4 


0.2-0.25 



kinetic (Ek) and magnetic (Em) energies: 0.9 < -E^ < 1-1, < Em < 0.05 for 
the first family, and 0.3 < Ek < 0.7, 0.07 < Em < 0.25 for the second one. 

The first family exists for < Rm < 58: the first attractor is the trivial 
steady state (5) (denoted by So) remaining stable up to^ Rm = 14. It becomes 
unstable in a pitchfork bifurcation, in which two stable mutually symmet- 
ric steady states (5*1^1 and 5*1^2) with a non-zero magnetic field emerge. The 
8-element symmetry group of Si^i is isomorphic to D4; however, it is distinct 
from Q. 

The steady states Si^i become unstable in a supercritical Hopf bifurcation 
at the interval 24 < Rm < 25. A stable periodic orbit Pi^i of period ~ 12 
appears in a vicinity of the steady state Si^i. The two orbits are mutually 
symmetric (they are interrelated hj h ). Each periodic orbit possesses all 8 
symmetries of its parent steady state, but each individual point of the orbit has 
a symmetry group of only 4 elements (it is isomorphic to D2). The orbits exist 
for 25 < Rm < 37 remaining attracting. A typical behaviour of a trajectory is 
shown on Fig. 2 for Rm = 30. After an initially small magnetic energy at first 
exponentially grows and afterwards decays with oscillations, it subsequently 
levels off (see a plateau about 400 time units long on Fig. 2a): the trajectory 
in the phase space evolves in a vicinity of Si^i (which is unstable now). In the 
further evolution the trajectory leaves the steady state, being attracted by Pi^i 
(see Fig. 2b). The orbits disappear in a subcritical Hopf bifurcation between 
Rm = 37 and Rm = 37.5 . For Rm = 37.5, 38 and 38.5 the two steady states 
Si^i are verified to be stable. 

Next bifurcation of Si^i is again a supercritical Hopf bifurcation. The 
emerging mutually symmetric periodic orbits P2i, i = 1,2 (one for each of 
the two steady states) possess all the symmetries of the steady states. For 
39 < Rm < 39.2 they are attracting, and their period is t{PI^) ^ 120. Be- 
havior of a sample trajectory in the phase space (see Fig. 3) resembles the 
one shown on Fig. 2; however, on Fig. 3a the plateau corresponding to the 
evolution in the vicinity of the unstable steady state Si^i virtually disappears, 
and the period is much larger. The orbit is located in a different region of the 
phase space (cf. Fig. 2b and Fig. 3b). At an Rm between 39.2 and 39.3 two 
mutually symmetric orbits P|j {i = 1, 2) of a twice larger period emerge (see 
Fig. 4a). A sequence of period-doubling bifurcations begins. The next period 
doubling occurs at the interval 39.4 < Rm < 39.45 (see period-four orbit P^j 
for Rm = 39.45 on Fig. 4b). At Rm = 39.5 a trajectory, initially close to Sq, 
is either already chaotic, or close to an orbit of a very long period P^^ (see 
Fig. 4c). This indicates that at Rm = 39.5 the period-doubling cascade is over. 
For Rm = 39.6 a trajectory with the same initial condition is attracted to a 
period-one orbit Pg^j. 

Another attractor exists for R„, = 39.6 . It reveals itself when the attractor 



^This is consistent with results of Galanti et al. (1992): since the ABC coefficients (6) 
are not normalised to satisfy A^ + B^ + C^ ~ ?, as they were in the cited paper, Reynolds 



numbers referred to here are lower by a factor (17/24) 
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for Rm = 39.7, which is unique in this region of the phase space (see below) 
is continued in smaller Rm (no attractor of this type is observed by further 
continuation in smaller Rm)- This is a periodic orbit P3, involving transitions 
between the two period-two orbits P21 and P22 (see Fig. 5a). The orbit is 
invariant under h and any symmetry from Q. In the course of temporal evo- 
lution, direction of magnetic field is reversed, behaviour of non-zero magnetic 
field Fourier coefficients is similar to that shown on Fig. 5b. A reversal takes a 
relatively short time, during which initially the energy of magnetic field attains 
the minimum and kinetic - the maximum, and subsequently magnetic energy 
blows up to a maximum, with the kinetic energy simultaneously reaching its 
minimum (Fig. 5c, d). Thereby, during a reversal magnetic energy is at first 
transformed into kinetic one, and then back into magnetic energy; for higher 
Reynolds numbers reversals exhibit such energy transformations as well. 

This bifurcation is apparently of the following nature. In the period- 
doubling cascade infinitely many unstable periodic orbits P.^^ {i = 1, 2) of 
periods t{P^^) = Mt[PI^) are created. Consider two such orbits P^ {m = 2 
in our case) related by h. Suppose that for some value of a control parame- 
ter the stable manifold of P™^ intersects the unstable manifold of P^, and a 
heteroclinic trajectory connecting these orbits emerges. Consider the symme- 
try, which maps one of the periodic orbits to another. Under this symmetry 
the heteroclinic trajectory is mapped to a distinct heteroclinic trajectory, and 
the two form a heteroclinic cycle. The cycle may give rise to a periodic orbit 
P3, when the control parameter varies. (This scenario is analogous to that 
in two-dimensional dynamical systems, where a homoclinic connection gener- 
ically implies for close values of the control parameter existence of periodic 
orbits in a vicinity of the homoclinic trajectory, and their period tends to in- 
finity when the critical value is approached.) Thus, the regime of reversals 
emerges via heteroclinic connection of two periodic orbits. 

For 39.7 < Rm < 41 sample trajectories, initially close to the trivial steady 
state 5*0, are attracted by a periodic orbit P4. While P3 (for R^ = 39.6) 
follows closely each of Pg^^ (Fig- 5a), this is not the case any more for P4 (for 
the higher Pm's; see Fig. 6). In agreement with our conjecture that emergence 
of the orbit is related to heteroclinic connection, the period decreases from 800 
for Rm = 39.7 to 500 for R^ = 41. 

For 39.6 < Rm < 41 the attractor involving reversals is unique. Behaviour 
of all magnetic field Fourier harmonics is coherent and resembles the one shown 
on Fig. 5b: time intervals, where a Fourier coefficient remains close to a cer- 
tain value, are separated by jumps changing the sign of the coefficient, but not 
its magnitude. In contrast with the coherent behaviour of all magnetic field 
Fourier harmonics for 39.6 < Rm < 41, they can be divided for 42 < Rm < 51 
into two groups with a qualitatively different behaviour. Behaviour of harmon- 
ics from the first group is the same as above, and hence the system exhibits 
reversals. Harmonics from the second group oscillate with a relatively short 
temporal period, some of their time averages do not vanish. The total time- 
averaged energy of harmonics from the first group is higher than that from 
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the second group (for instance, cf. Fig. 7a and 7b). The Earth's magnetic 
field shows a similar behavior: Spherical harmonics can be categorized into 
the dipole and the quadrupole families according to their symmetry about the 
equatorial plane (Roberts and Stix, 1972), and paleomagnetic data suggest 
that a magnetic reversal in the Earth is more likely to occur, when the ratio 
of the magnitude of the quadrupole-family component of the magnetic field 
to that of the dipole-family component is high (Merrill and McFadden, 1988). 
For 42 < Rm < 51 our attractors have no the symmetry h on average (since 
the time average of some harmonics is non-zero), and thus h maps an attractor 
to a different attractor with reversals. Again, a parallel can be drawn with a 
property of the magnetic field of the Earth: Paleomagnetic data indicate that 
the time-averaged magnetic field of the Earth lacks the north-south symmetry 
and there are significant (though small) differences between the states of the 
normal and reverse polarity (Merrill et ai, 1979), this implying the absence of 
the symmetry h on average. 

For Rm = 42 and 43 the regime is quasi-periodic: the orbit P4 bifurcates 
into two attracting tori Ti^j {i = 1,2) interrelated by h (see Fig. 7). The 
second frequency can be observed on Fig. 7a, showing real part of the Fourier 
component bQnoi magnetic field. The time average of this Fourier component 
is negative (Fig. 7a); hence, the symmetry h maps this attractor to a distinct 
one. (For 37.5 < Rm < 41 the component vanishes, because former attractors 
[P4 for 39.7 < i?„ < 41, P3 for R^ = 39.6, P2™ for 39 < P„ < 39.6 and Si^i 
for 37.5 < Rm < 38.5] possess the symmetry hs^ = hsiS2-) 

For Rm = 45 two mutually symmetric attractors persist (Fig. 8): each 
torus Tii bifurcates into a chaotic attractor C2,j. Reversals are less regular, 
yet they are too ordered in comparison with those of the Earth's magnetic 
field. 

For Rm = 46 (Fig. 9) behaviour of a sample trajectory suggests existence 
of four (unstable) steady states 82,1 (see plateaux at 1500 < t < 2100). They 
are mutually symmetric, and have the same 4-element symmetry group gen- 
erated by hsi and S3. Apparently two 6*2,1 's emerged from each Si^i in a 
pitchfork bifurcation. Initially (at < t < 700) the sample trajectory for 
Rm = 46 undergoes several reversals, similar in nature to those observed for 
39.7 < Rm < 45. Afterwards it is attracted to a steady state S2,i. The be- 
haviour near the steady state resembles the Shilnikov attractor: the trajectory 
approaches the steady state along the two-dimensional eigenspace associated 
with the complex eigenvalues with the maximal negative real part, and leaves 
it along the one-dimensional unstable manifold. Subsequent sample evolution 
consists of transitions between the steady state 6*2,1 and the region of the phase 
space where a former chaotic attractor C2,j was located (duration of repeat- 
ing events in the saturated regime is ~ 1700). In particular, large-amplitude 
excursions at 2350 < t < 2500, 2300 < t < 2450, 4000 < t < 4150, etc., are 
reminiscent of the behaviour of a C2,j trajectory for Rm = 45. Thus, appar- 
ently two former chaotic attractors C2,i have disappeared in collision with the 
four steady states 5*2, j to give rise to four new chaotic attractors Ca^j. It is 
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notable that the sample trajectory leaves a vicinity of the 82,1 in alternating 
directions along the unstable manifold. For Rm = 47 the behaviour is similar 
to the one at Rm = 46. 

For Rm = 48 the system possesses four mutually symmetric attractors 
C^^i. Events are longer than in the previous regime, about 3000 time units 
(see Fig. 10). They exhibit a new feature, a phase of initially exponentially 
decaying oscillations (e.g. at 1700 < t < 2300 and 4500 < t < 5200). This 
indicates existence of a weakly unstable periodic orbit of a new kind, P5 j 
{i = 1, 4), in a vicinity of the former chaotic attractor C^^i. A similar behaviour 
is also observed for Rm = 49. 

For Rm = 50, 51 (Fig. 11) the four C4^j are superceded by two new attrac- 
tors, Cs^j. Large amplitude oscillations, e.g. at 7500 < t < 8500, are due to 
attraction of the trajectory by weakly unstable tori T2^j {i = 1,4), which have 
bifurcated from the periodic orbits P^^i. In the saturated regime behaviour 
consists of three phases: (i) a trajectory is close to a periodic orbit P^^i (e.g. 
at 2300 < t < 3000 and 4300 < t < 5000); (ii) the trajectory evolves in the 
vicinity of a torus T2,i (e.g. at 3200 < t < 4000 and 5300 < t < 6000); {in) the 
trajectory abruptly jumps toward the second periodic orbit P^^i' to reproduce 
the sequence of phases. 

For Rm = 52 each of the two former chaotic attractors splits into two 
attracting tori T2^i (Fig. 12; cf. Figs, lib and 12b). The four new attractors 
are mutually symmetric, they are stable for 52 < Rm < 58. 

For 57 < Rm < 60 a sample trajectory initially close to 5*0 is attracted 
by the torus T3, belonging to the second family. The family emerges between 
Rm = 15 and Rm = 16: in addition to the steady states Si^i, for Rm > 16 
the system possesses another attractor, a stable periodic orbit Pq. The orbit 
does not exist for Rm < 15; apparently it appears in a saddlenode bifurcation. 
For Rm = 26 the orbit is unstable and a torus T3 has appeared in a Hopf 
bifurcation. T3 remains an attractor for all higher considered Rm (specifically, 
for Rm = 26, 27, 30, 40, 57 and 60). For Rm > 57 the torus attracts trajecto- 
ries (in the phase space) initially close to the trivial steady state So, and for 
Rm > 59 it is the only global attractor of the system. The torus T3 and the 
orbit Pq were traced back from Rm = 57 by continuation in Rm- 

Studying numerically a dynamical system, one can never guarantee that 
all attractors are found. The computations are intensive (a run with the reso- 
lution 32^^ Fourier harmonics up to t = 5000 has required several dozens hours 
of one processor), prohibiting to carry out more detailed computations. For 
instance, if attractors of two or more new types emerge for a particular con- 
sidered Rm, in principle we could have missed some of them; however, it is 
unlikely that several attractors of different morphology are born at the same 
or close values of the magnetic Reynolds number. In our estimations all at- 
tractors, whose basin of attraction contains the trivial steady state, have been 
identified - up to a natural limitation: since a period-doubling cascade involves 
an infinite number of bifurcations, it is impossible to capture numerically all 
bifurcations within the cascade. Attractors, which do not attract small per- 
il 



turbations of the trivial steady state and which are far from it in the phase 
space, (if such attractors exist in our system) are not involved in the scenario 
of emergence of magnetic field reversals, and they are out of the scope of the 
present investigation. 



5 Attractors of the MHD system for R^ = 40 

Simulations are performed for Rm = 40 and for several values of R with 
an initial condition, same for all runs, being a small perturbation of 5*0 (see 
Table 2). The question we address is how the temporal behaviour changes with 
R, in particular, if for higher R reversals persist, or whether magnetic field 
vanishes in a saturated regime. We do not attempt to identify all bifurcarions 
occurring when R varies. Computations with the same R and the same initial 
condition for the flow are also performed for the purely hydrodynamic system 
(2. a), in order to compare the behaviour and attractors of the two systems. 
(If (2) possesses attractors with b = different from the trivial steady state, 
an initially growing magnetic field can die out in the nonlinear regime, the 
resulting attractor being purely hydrodynamical, see Brummell et al, 1998; 
Matthews, 2000.) 

For R = 3 the detected attractor is the torus Ti (which is an attractor for 
R = 4 and 26 < Rm < 60). For R = 6 reversals take place, similar to those 
observed for R = 4 and Rm = 45, and the attractor is C2,i. Thus, for R close 
to 4 no new attractors are found. 

For i? = 10 a new attractor was found in the full MHD system - a peri- 
odic orbit Pj with a symmetry group of 8 elements. Comparison of Tables 2 
and 3 reveals no relation of attractors of the hydrodynamic and MHD systems 
for R = 10, 15 and 20. For these Reynolds numbers a small magnetic field 
(magnetic energy Em is below 0.07; see Table 2) drastically changes behaviour 
of the system - the hydrodynamic and MHD systems have attractors of dif- 
ferent types, and the average kinetic energy Ek decreases significantly, e.g. 
Em = 0.01 and E^. decreases from 0.9 (the hydrodynamic case) to 0.7 (the 
MHD case) for R = 20. 

For i? = 15 behaviour of the sample trajectory of the MHD system in 
saturated regime is chaotic. But before the saturated regime sets in, the 
trajectory is attracted by Pj and for 1000 < t < 2000 remains close to this 
periodic orbit (see Fig. 13), which is now weakly unstable. For R = 20 the 
trajectory with the same initial condition is attracted to a new torus T4 j; 
initially the temporal behaviour of the trajectory resembles the one observed 
for R = 15 (cf. Figs. 14 and 14). 

For R = 25 the MHD system possesses new chaotic attractors Cj^i {i = 1,2) 
with a trivial symmetry group. They resemble the chaotic attractor Cq ob- 
served for i? = 15 (cf. Fig. 13 and Fig. 15), but unlike Cj^i, Cq is unique and 
has a symmetry group of four elements. Comparison of Cj^i with the attractor 
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Table 2. Attractors, detected for the MHD system with the force (3), (4), (6) 
for Rm = 40 and 3 < i? < 25. The third column shows the number of elements 
of the symmetry group for which an attractor is pointwise invariant, and the 
fourth - generators of the group. 



R 


Attractors 


Number of 
symmetries 


Generators 


Ek 


Em 


R = 3 


n 


4 


hS2 


0.4 


0.25 


R = A 


periodic orbit P4 


8 


hsi, S2 


0.95 


0.03 


R = 6 


chaotic C2,i {i = 1, 2) 


4 


hsi, S3 


0.95 


0.03 


R=10 


periodic orbit Pj 


8 


hsi, S2 


0.5 


0.06 


R=15 


chaotic Ce 


4 


Sl, S3 


0.6 


0.04 


R = 20 


tori T4,, {t = 1, 2) 


4 


Sl, S3 


0.7 


0.01 


R = 25 


chaotic C-i^i {i = 1, 2) 


1 


e 


0.5 


0.02 



Table 3. Attractors, detected for the hydrodynamic system (2. a) with the force 
(3), (4), (6) for 3 < i? < 25. The third column shows the number of elements 
of the symmetry group for which an attractor is pointwise invariant, and the 
fourth column - generators of the group. 



R 


Attractors 


Number of 
symmetries 


Generators 


Ek 


i? = 3,4,6 


uabc 


8 


Sl, S2 


1.06 


R= 10* 


'^ABC 


8 
8 


Sl, S2 
Sl, S2 


1.06 
0.99 


i?= 15 


S^ 


8 


Sl, S2 


0.93 


R = 20 


chaotic C'l 


4 


S2 


0.9 


R = 25 


chaotic C2 


1 


e 


0.7 



* Sl for i? = 10 was traced back from R = 15; it is not observed for R = 6. 
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of the non-magnetic Navier-Stokes equation, which also is chaotic, reveals that 
the influence of magnetic field is in some sense stabilising: Fourier coefficients 
of the flow experience fewer jumps and the amplitude of their oscillations is 
much smaller (compare Fourier components of flows on Figs. 15b and 16b). An 
exponential growth of the initially small magnetic field in the sample evolution 
of the system begins only at t = 200 (Fig. 15c). Accordingly, until t = 300 the 
flow evolution is similar to the one in the absence of magnetic field. Magnetic 
field starts growing, when the trajectory is attracted by a weakly unstable pe- 
riodic orbit in the hydrodynamic subspace. Departure from this orbit causes 
an initial decay of magnetic field before the onset of saturated behaviour. 

For R = 15 and 20 behaviour of Fourier components of magnetic field is 
similar to the one shown on Fig. 15d. It consists of chaotic irregular small- 
period oscillations about zero, different from the behaviour observed for R = 4. 
Thus for Rm = 40, 15 < i? < 25 the considered system is not in a regime 
resembling reversals. 

6 Conclusion 

The bifurcation scenario leading to emergence of reversals, which we put for- 
ward, proves feasible. Reversals are found for small R, for which the hydro- 
dynamic system has a unique globally stable steady state. Complexity of the 
sequence of bifurcations obtained in simulations is comparable to that of the 
hydrodynamic system examined by Podvigina and Pouquet (1994) and Pod- 
vigina (1999). Numerous types of behaviour are identified, which involve a 
large variety of attractors and unstable invariant sets of various types affect- 
ing behaviour of evolutionary solutions. The regime of reversals arises as a 
result of merging of two distinct attractors, apparently via heteroclinic con- 
nection due to intersection of stable and unstable manifolds of periodic orbits 
(in contrast to the system of Armbruster et al. (2001), where a heteroclinic 
cycle emerged simultaneously with magnetic steady states). 

Other regimes with a non-zero magnetic field are observed when any of the 
Reynolds numbers is increased. Thus in our computations reversals are not so 
robust, as in simulations of Glatzmaier and Roberts. The range of parameters 
considered here is rich in bifurcations: for i? = 4 we found 14 bifurcations in 
the interval 14 < Rm < 60 (counting the infinite number of period-doubling 
bifurcations as one), 7 distinct types of attractors involve magnetic field rever- 
sals. In another dissimilarity with the results of Glatzmaier and Roberts (and 
with the Earth's magnetic field) in our simulations Em < E^ (cf. E^/ E^ ~ 10^ 
for the geodynamo). In our computations reversals appear more regular than 
those in the Earth, since much smaller values of Reynolds numbers are con- 
sidered. The dissimilarity with reversals in the Earth also owes to the fact 
that we solve a simplified system of PDE's (neither an account of thermal or 
sedimentation-driven convection is taken, nor that of the Coriolis force) in an 
idealised space-periodic geometry. Nevertheless, it appears possible to mim- 
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ick the presence of two time scales (average periods of constant polarity are 
much longer than the average time of a reversal itself); instances of excursions 
of magnetic field (i.e., incomplete reversals) can be observed as well (see e.g. 
Figs. 8a and 10). 

Our results show that magnetic field reversals are an inherent feature of 
a nonlinear dynamical system of the MHD type. The fundamental reason for 
their existence is the symmetry /i, changing polarity of the magnetic field. 
Neither external triggering random processes, such as the ones postulated in 
earlier studies, nor symmetry breaking due to inhomogeneity of the outer core 
boundaries, nor rotation, nor a specific heteroclinic structure or symmetries of 
the flow sustaining a non-linear dynamo are required to give rise to a sequence 
of reversals. Such physical processes, as heat transfer or sedimentation are 
also unnecessary per se (as soon as the system is under the action of any body 
force providing input of energy) . 
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Figure 1. Bifurcation diagram of the detected MHD attractors for R = 4 and 
< Rm < 58. For each type of attractors the interval of -R^'s (horizontal axis, 
non-uniform scale) is indicated, for which attractors of this type have been 
detected in computations. Labeling of attractors is explained in Section 4 (see 
also Table 1). Only one representative symmetry-related attractors is shown. 
The attractor P21 for Rm = 39.6 and attractors from the second family are 
not shown. Stable steady states are represented by solid lines, the unstable 
ones - by dashed hues. 
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Figure 2a. Magnetic energy (vertical axis) as a function of time (horizontal 



axis) for i? = 4 and Rm = 30. 
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Figure 2b. Projection of the trajectory in saturated regime (a periodic orbit 
Pi^i) on the plane of Fourier coefficients Imfeg ^ 2 (horizontal axis) and RcWq ^ 2 
(vertical axis) for R = 4 and Rm = 30 (same run as on Fig. la). 
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Figure 3a. Magnetic energy (vertical axis) as a function of time (horizontal 



axis) for i? = 4 and Rm = 39.2. 
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Figure 3b. Projection of the trajectory in saturated regime (a periodic orbit 
P2i) on the plane of Fourier coefficients Imfeg ^ 2 (horizontal axis) and RcWq ^ 2 
(vertical axis) for R = 4 and Rm = 39.2 (same run as on Fig. 2a). 
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Figure 4a. Projection of the trajectory in saturated regime (a periodic orbit 
P2i) on the plane of Fourier coefficients ImfeQ ^ 2 (horizontal axis) and Re^Q ^^ 2 
(vertical axis) for R = 4 and Rm = 39.4. 
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Figure 4b. Projection of the trajectory in saturated regime (a periodic orbit 
P2i) on the plane of Fourier coefficients Imfeg ^ 2 (horizontal axis) and RefQ ^^ 2 
(vertical axis) for R = 4 and Rm = 39.45. 
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Figure 4c. Projection of the trajectory in the phase space in saturated regime 
(a periodic orbit P^^) on the plane of Fourier coefficients ImfoQ ^ 2 (horizontal 
axis) and Refg ^^ 2 (vertical axis) for R = 4 and Rm = 39.5. 
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Figure 5a. Projection of the trajectory in saturated regime (periodic orbit P3) 
on the plane of Fourier coefficients Im 6j 1 2 (horizontal axis) and Re f g 1 2 (ver- 
tical axis) for R = 4 and Rm = 39.6. 
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Figure 5b. Fourier coefficient Im6Q^2 (vertical axis) as a function of time 
(horizontal axis) for R = 4 and Rm = 39.6 (same run as on Fig. 4a). 
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Figure 5c. Magnetic energy (vertical axis) versus time (horizontal axis) for 
R = 4 and Rm = 39.6 (same run as on Fig. 4a). 
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Figure 5d. Kinetic energy (vertical axis) versus time (horizontal axis) for 
R = 4 and Rm = 39.6 (same run as on Fig. 4a). 
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Figure 6. Projection of the trajectory in saturated regime (periodic orbit P4) 
on the plane of Fourier coefficients Im 6j 1 2 (horizontal axis) and Re Wq ^ 2 (ver- 
tical axis) for R = 4 and Rm = 39.7. 



25 



0.010 



0.005 



D.OOO - 



-0.005 



-0.010 




3000 



3200 



3400 



3600 



3800 



4000 



Figure 7a. Fourier coefficient Im &J i 2 (vertical axis) as a function of time 
(horizontal axis) for R = 4 and Rm = 42. 
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Figure 7b. Fourier coefficient Re 6q ^^ ^ (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 42 (same run as on Fig. 6a). 
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Figure 8a. Fourier coefficient Im &J i 2 (vertical axis) as a function of time 
(horizontal axis) for R = 4 and Rm = 45. 



0.015 



.010 - 



0.005 - 



0.000 irjlk 



-0.005 - 



-0.010 - 



-0.015 



-0.020 L 




I I'V 



ft.l1llll''llliil,l'jt,'i,ll,l 



1000 



2000 



3000 



Figure 8b. Fourier coefficient Re 60 ^^ ^ (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 45 (same run as on Fig. 7a). 
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Figure 9a. Fourier coefficient Im &J i 2 (vertical axis) as a function of time 
(horizontal axis) for R = 4 and Rm = 46. 
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Figure 9b. Fourier coefficient Re 6q ^^ ^ (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 46 (same run as on Fig. 8a). 
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Figure 9c. Fourier coefficient Re 6^ o i (vertical axis) as a function of time 
(horizontal axis) for R = 4 and Rm = 46 (same run as on Fig. 8a). 
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Figure 10. Fourier coefficient Im6j]^2 (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 48. 
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Figure 11a. Fourier coefficients Imfeg^g (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 51. 
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Figure lib. Projection of the trajectory in saturated regime on the plane of 
Fourier coefficients Im 6q ^ 2 (horizontal axis) and Re f g ^ 2 (vertical axis) for 
R = 4 and Rm = 51 (same run as on Fig. 10a). 
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Figure 12a. Fourier coefficient lm.bl^2 (vertical axis) as a function of time 
(horizontal axis) for i? = 4 and Rm = 52. 
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Figure 12b. Projection of the trajectory in saturated regime (torus T^^i) on 
the plane of Fourier coefficients Im 6q ^ 2 (horizontal axis) and Re f g ^ 2 (vertical 
axis) for R = 4 and Rm = 52 (same run as on Fig. 11a). 
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Figure 13. Magnetic energy (vertical axis) as a function of time (horizontal 
axis) for R = 15 and Rm = 40. 
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Figure 14. Magnetic energy (vertical axis) as a function of time (horizontal 
axis) for R = 20 and Rm =40. 
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Figure 15a. Kinetic energy (vertical axis) as a function of time (horizontal 



axis) for R = 25 and Rm = 40. 
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Figure 15b. Fourier coefficients of the flow (vertical axis: solid line - Ref^oo; 
dashed line - Refg^Q, dot line - ReWgQ^) as a function of time (horizontal 
axis) for R = 25 and Rm = 40 (same run as on Fig. 14a). 
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Figure 15c. Magnetic energy (vertical axis) as a function of time (horizontal 
axis) for R = 25 and Rm = 40 (same run as on Fig. 14a). 
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Figure 15d. Fourier coefficient Im6Q^2 (vertical axis) as a function of time 
(horizontal axis) for R = 25 and Rm = 40 (same run as on Fig. 14a). 
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Figure 16a. Kinetic energy (vertical axis) as a function of time (horizontal 
axis) for i? = 25 in the hydrodynamic problem (the initial condition for the 
flow is the same as on Fig. 14). 
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Figure 16b. Fourier coefficients of the flow (vertical axis: solid line - Reff qq, 
dashed line - Refg^g, dot line - Re^gg^) as a function of time (horizontal 
axis) for R = 25 (same run as on Fig. 15a). 
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